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ABSTRACT 

We perform a maximum likelihood analysis of the cluster abundance measured in the SDSS using 
the maxBCG cluster finding algorithm. Our analysis is aimed at constraining the power spectrum 
normalization er 8 , and assumes flat cosmologies with a scale invariant spectrum, massless neutrinos, 
and CMB and supernova priors £l m h 2 = 0.128 ± 0.01 and h = .72 ± 0.05 respectively. Following 
the method described in the companion paper IRozo et all (|2007f ). we derive a% = 0.92 ± 0.10 (la) 
after marginalizing over all major systematic uncertainties. We place strong lower limits on the 
normalization, erg > 0.76 (95% CL) (> 0.68 at 99% CL). We also find that our analysis favors 
relatively low values for the slope of the Halo Occupation Distribution (HOD), a = 0.83 ± 0.06. 
The uncertainties of these determinations will substantially improve upon completion of an ongoing 
campaign to estimate dynamical, weak lensing, and X-ray cluster masses in the SDSS maxBCG cluster 
sample. 

Subject headings: cosmology: theory — cosmological parameters — galaxies: clusters — galaxies: 
halos 



in 
m 
O 

o 

6 
c3 



1. INTRODUCTION 

One of the most important problems in observational 
cosmology today is to resolve the question of whether or 
not dark energy takes the form of a cosmological con- 
stant. While current geometric probes of dark energy 
such as supernovae and baryon acoustic oscillations un- 
equivocally tell us that dark energy exists, a complemen- 
tary probe of the dark energy evolution can help us dis- 
tinguish between a cosmological constant and a dynam- 
ical da rk energy. The growth of structure is one such 
prob e (lEke et al.lll996t iHolder et al.l l2001l: lEvrard et all 
120021 iMolnar et al.ll2004f ). In particular, given a geomet- 
ric probe, an accurate determination of the amplitude 
of the power spectrum at two different times directly 
constrains the growth between the two epochs, and can 
in principle help determine if the dark energy density 
evolves with redshift. We know the power spectrum am- 
plitude at the time of last scattering with high preci- 
sion thanks to Cos mic Microwave Back ground (CMB) 
experiments (see e.g. ISpergel et al.ll2006l and references 
therein), so a precise determination of the current power 
spectrum normalization may, in principle, distinguish be- 
tween a dynamical dark energy component and a cosmo- 
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logical constant. 

It is well known that the abundance of massive ha- 
los in the local universe depends strongly on the ampli- 
tude of the matter power spect rum as- 11 In particular , 
from theoretical considerations (Press fc Schechterlll974l : 
iBond et al.lll99lUSheth fc Tormenll2002fl one expects the 
number of massive clusters in the universe to be exponen- 
tially sensitive to ag, a picture that has been confirmed 
with extensive numerical simulations. Consequently, the 
number of galaxy clusters within a given survey region 
ought to be able to provide powerful constraints on as, 
and, indeed, the literature is rife with these type of stud- 
ies. 

Unfortunately, cluster abundance determinations of as 
have to overcome a large variety of difficulties. For in- 
stance, the abundance of massive halos in the universe is 
sensitive not only to as but also to Q m , the mean matter 
density of the universe, which implies that constraints 
on the number density of massive halos typically result 
in large degener acies between fl m and as (though see 
IRozo et al.ll2004j ). In fact, the main the main obstacle to 
accurate as measurements is syste matic uncertainties in 
mass estima tes of clusters (see e.g. iPierpaoli et aill2003l : 
iHenrvl l2004f ). Consequently, new analyses that prop- 
erly marginalize over such systematic uncertainties are 
of particular importance to interpret cluster abundance 
constraints within a broad cosmological context. 

In this work, we use the techniques developed in 
I Rozo et "ail (|2007t ) to analyze the Sloan Di gital Sky Sur- 
vey (SD SS) maxBCG cluster sample from lKoester et aTl 
(|2007aj ) . The analysis of these data uses information 
from the lRozo et aLl(|2007h companion paper on the form 
of the maxBCG selection function, which connects our 
observable mass proxy — the cluster richness — to halo 
mass. As discussed in in that work, at this time our 
understanding of the maxBCG selection function is in- 

11 Here, we characterize the present day amplitude of the power 
spectrum with the usual parameter as, the rms amplitude of den- 
sity perturbations in spheres of 8/i —1 Mpc radii. 
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complete, which implies that we have large systematic 
uncertainties in the selection function. Nevertheless, the 
large mass range probed by the maxBCG cluster sam- 
ple allows us to marginalize over these uncertainties and 
still recover competitive estimates for erg, albeit with the 
inclusion of cosmological priors for fl m h 2 and h in our 
analysis. Importantly, the results presented in this work 
ought to be interpreted as an upper limit of how well 
we can expect erg to be constrained from the current 
maxBCG cluster sample in the near future. Not only 
will we soon be able to include additional data such as 
weak lensing and dynamical cluster mass estimates, but 
we also expect our understanding of the maxBCG cluster 
selection function to improve both as an expanded suite 
of simulations used to calibrate the maxBCG selection 
function become an even more accurate approximation 
to reality and as we work towards a more robust richness 
estimator. 

The layout of the paper is as follows. In [21 we briefly 
describe the maxBCG cluster finding algorithm used to 
identify the cluster sample analyze d in this work. jp3 
summarizes the model described in iRozo et ail (|2007l )T 
including a description of the various parameters and 
priors used in our analysis. Our cosmological constraints 
are presented in and discussed in detail in !j5] 

2. THE MAXBCG CLUSTER-FINDING ALGORITHM 

Details of how the maxBCG cluste r -finding algorithm 
works can be found in iKoester et al.1 (|2007al) . Here, we 
only summarize the main elements of the cluster-finding 
algorithm. 

MaxBCG is an optical cluster-finding algorithm that 
relies on photometric measurements to overcome projec- 
tion effects. To detect clusters, maxBCG uses the well 
known observational fact that galaxy clusters contain 
a large number of so called ridgeline galaxies: bright, 
red, early type galaxies that populate a narrow ridge- 
line in color-magnitude space as a function of redshift. 
The color distribution of these galaxies is modeled as 
a narrow Gaussian, while their two dimensional spa- 
tial distribution about the cluster center is modeled as 
a projected Navarro , Frcnk, and White (NFW) profile 
(| Navarro et alj fl997l ) . To determine the cluster center, 
maxBCG relies on the observational fact that, in the 
vast majority of clusters, there is a clear Brightest Clus- 
ter Galaxy (BCG) at or near the center of the cluster. 
These BCG galaxies tend to be extremely luminous and 
red, populating the tip of the color-magnitude ridgeline. 
Their luminosity and colors are also modeled as narrow 
Gaussians. 

Given our model, one can then compute the likelihood 
that a particular galaxy is the BCG galaxy of a cluster 
by computing the product Cbcg^-r, where Cbcg{ z ) is 
the likelihood that the galaxy under consideration has 
the observed colors and magnitude assuming that it is 
a BCG at redshift z, and Cr is the likelihood that the 
galaxy distribution around the the candidate BCG will 
occur under the assumption that the BCG is at the cen- 
ter of cluster (though this is unlikely t o be universally 
true; see e.g. Ivan den Bosch et~aT]|2005D . In computing 
the ridgeline likelihood Cr, only galaxies found within 
3 Mpc of the candidate BCG and in a narrow (3c) color 
window around the expected color for ridgeline galax- 
ies at that redshift are considered. The total likelihood 



C = Cbcg^-r is evaluated at a grid of rcdshifts, and a 
photometric redshift estimate for the cluster is obtained 
by maximizing the likelihood function. 

The end result of this process is a likelihood assign- 
ment for every galaxy in the survey. The galaxy list is 
then rank ordered according to likelihood, and the most 
likely BCG is selected as a BCG. A first rough richness 
measurement is made by counting the number of galax- 
ies within 1 Mpc, which is then used to estimate a char- 
acteristic rad i us for the cluster using the results from 
lHansen et al.1 (|2005l ) . All ridgeline galaxies above a lumi- 
nosity cut of 0.4L* within this scale radius are considered 
cluster members, and the number of cluster members is 
define d as the cluster's rich ness, denoted here by N a b s 
and in IKoester et~aTI (|2007aD by N'™° 12 . All BCG can- 
didates within a galaxy overdensity of 200 and within 
Az = 0.02 of the cluster just found are dropped from 
the candidate BCG list. The procedure is then iterated, 
resulting in a cluster catalog where each cluster has an 
assigned cluster center, a members' list, a richness, and 
a photometric redshift estimate. 

3. THE MODEL 

In this work, our main observable is the number of 
clusters within the survey volume as a function of cluster 
richness, i.e. the richness function. We now summarize 
the basic picture behind our analysis. A d etailed pre- 
sentat ion of this formalism can be found in IRozo et alJ 

(Hqq3). 

3.1. The Model at a Glance 

Suppose we wish an expression for the number of clus- 
ters of a given richness. In general, if P(N i, s \m) is the 
probability that a halo of mass m is detected as a cluster 
with N f,s galaxies, the mean density of these clusters is 
simply 

( n ) = [dm ^P-P(N obs \m). (1) 
J am 

The main idea behind our analysis is to split P(N b s \m) 
into two: there is a probability P(Nx\m), namely the 
Halo Occupation Distribution or HOD, that determines 
the intrinsic scatter between a halo's mass and its rich- 
ness, and a second distribution P(N i, s \Nt) that charac- 
terizes the observational scatter. If c(Nt) is the proba- 
bility that a halo with Nt galaxies is detected, then the 
probability that a halo of mass m is detected as a cluster 
with N bs galaxies is 

P(N ob s\m) = Y,<N T )P{N ob s\N T )P(NT\m). (2) 

Note that, by definition, c(Nt) is also the expected frac- 
tion of detected halos, and hence we refer to it as the 
completeness. Finally, suppose p(N b s ) is the probabil- 
ity of a cluster with N b s galaxies being a real detection, 
i.e. of the cluster corresponding to an actual halo. By 
definition, p(N b s ) is also the expected fraction of clus- 
ters that are real, and hence we call p(N b s ) the purity 
of the sample. If p(N b s ) ^ 1, then the observed number 
of clusters will be boosted relative to our above estimate 

12 For technical reasons, in this work the luminosity cut for Nt 
is defined in the r-band, whereas N a i, a has an i-band luminosity 
cut. 
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Fig. 1. — Effec tive mass select i on fun ction for the maxBCG SDSS 
cluster sample of Kocstcr ct al. (2007a) in our maximum likelihood 
model. The thin solid lines are the selection function for each of 
our 10 richness bins. The selection function for the richest bin 
continues to grow at larger ma sses due to the ex trapolation of the 
decreasing purity functionfsce IRozo et al]|2007l) . Note, however, 
that all abundances remain finite since the halo mass function is 
decreasing exponentially fast at these mass scales. The thick, solid 
line is the net sum of the individual bins, and gives the effective 
mass selection for the cluster sample as a whole. 

by a factor l/p(N b s ), so the mean number density of the 
clusters becomes 

(n) = [dm ^P-i>(m) (3) 
J dm 

where 

= ^r^J2 c ( N ^ P ( N ^WT)P(N T \m). (4) 

The quantity ip(m) is the effective mass binning of our 
cluster sample, that is, ip(m) is the fraction of mass m 
clusters that fall into the richness bin of interest. Note, 
however, that when p(N b s ) < 1, then ip(m) can be larger 
than unity. The true mass binning is obtained by setting 
p(N bs) = 1- Figure Q] shows the effective mass binning of 
of our cluster sample assuming the maximum likelihood 
values of all relevant model parameters as determined by 
our analysis (see below). 

3.2. The Likelihood Function 

The above model for cluster abundances tells us what 
the expected cluster abundance will be. For our analy- 
sis, however, we wish to know the likelihood of observ- 
ing a particular data set given s ome cosmological and 
HOD parameters. As described in lRozo et al.l (|2007( ). we 
model the probability of observing a realization given a 
set of cosmological and HOD parameters as a Gaussian. 
While more accu rate likelihood functions can be found 
in the literature (|Hu fc Cohnll2006t lHolderll20060 , these 
ignore correlations due to scatter in the mass-obscrvablc 
relation, and thus we have opted for a simple Gaussian 
model, which is expected to hold if bins are sufficiently 
wide to include a large (> 10) number of clusters. 

The contributions to the correlation matrix that we 
consider are: 

1. A Poisson contribution due to the Poisson fluctua- 
tion in the number of halos of mass m within any 
given volume. 



2. A sample variance contribution reflecting the fact 
that the survey volume may be slightly overdense 
or underdense with respect to the universe at large. 

3. The bin-to-bin scatter arising from the stochastic- 
ity of N b s as a function of m. 

4. A contribution due to the statistical uncertainties 
associated with photometric redshift estimation. 

5. A contribution due to the stochastic nature of the 
completeness and purity functions. That is, if we 
know the expected purity and completeness of an 
infinite sample, any finite sample may have slightly 
different fractions of true detections and false pos- 
itives. 

The detailed construction of th is likelihood function 
can be found in IRozo et ail (|2007f ). For the purposes of 
this work, the most important aspe ct of our likelihoo d 
function is that, as demonstrated in IRozo et al.l (|2007f) . 
our likelihood analysis correctly recovers the underly- 
ing cosmological and halo occupation distribution pa- 
rameters to within the intrinsic degeneracies of the data 
for mock maxBCG catalogs constructed using the tech- 
niques in (Wechslcr ct al. 2007, in preparation). Con- 
sequently, we are fully confident that our analysis tech- 
nique is sound, robust, and that it properly takes into 
account the various sys tematic uncertainties that affect 
the construction of the iKoester et al.l fct )()";>;■ maxBCG 
catalog. 

3.3. Model Parameters 

For reference, we summarize here all of the relevant 
parameters for our model. The cosmological parame- 
ters we considered are <78,Q m , and h. The power spec- 
trum is taken to be scale invari ant, and we use t he low 
baryon transfer function from lEisenstein fc Hul (|1999f ) 
with zero neutrino masse s. The baryon densit y is fixed 
to the WMAP3 value in iSpergel et all (120061). All re - 
sults reported in this work use a IJenkins etaTl (|200lh 
mass function, and we have explicitly checked that the 
expectation values for all parameters of interest recov- 
ered using differen t mass function parameterizations (in 
parti cular those of ISheth fc Tormenl 120021 : 1 Warren et al.l 
|2005[ ) f all well within the la error bars recovered us- 
ing the IJenkins et al.l (|2001| ) mass functi on. We also as- 
sume a flat AC DM universe 13 . FollowinglKraytsov et al.l 
(|2004L see also lZheng eraTll2005t lYang et al.ll2005al ) we 
assume that the total number of galaxies in a halo takes 
the form Nt = 1 + N sat where N sat , the number of satel- 
lite galaxies in the cluster, is Poisson distributed at each 
halo mass m with an expectation value (N sat \m) given 

by 

(N sat \m)= C^-Y- (5) 

13 Note that because of our cluster abundance determination of 
<T8 is both local and uses only a narrow redshift range, the con- 
straints we recover from the sample are largely independent of the 
dynamics of dark energy. Consequently, we do not expect assum- 
ing a ACDM universe will bias our results in any significant way. 
As error bars shrink, however, this assumption will undoubtedly 
need to be relaxed. 
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Here, M\ is the characteristic mass at which halos ac- 
quire one satellite galaxy. Note that in cluster abun- 
dance studies, the typical mass scale probed is consider- 
ably larger than M\. Nevertheless, the above parameter- 
ization is convenient because degeneracies between HOD 
and cosmological parameters take on particu larly simple 
forms when parameterized in this way (see iRozo et al.l 
12001 . 

Our model also includes a large number of nuisance 
param eters. These are described in detail in lRozo et al.l 
12001 . There, we demonstrate that the completeness 
function is flat as a function of richness, and hence can 
be described with only one nuisance parameter c. The 
purity function, on the other hand, is clearly peaked, 
with purity decreasing both in the high and low richness 
limits. We found that we can accurately describe the 
purity function with two nuisance parameters, pa and 
Pi. Two more parameters, Bq and /3, describe the mean 
value of N i, s for clusters at fixed Nt- These parameters 
characterize the amplitude and slope of the mean rela- 
tion {N i,s\Nt) respectively. An additional parameter B\ 
describes the variance of N b s at fixed Nt, which is taken 
to be simply proportional to (A,,^! Ax) 14 . The form of 
the distribution P(N i, s \Nt) is taken to be a discrctized 
Gaussian based on simulations (see below). Finally, two 
additional nuisance parameters (b z ) and o~ z calibrate the 
bias and scatter of our photometric rcdshift estimates. A 
summary of the equations that define all of our nuisance 
parameters is presented in Appendix [AJ 

3.4. Model Priors 

In IRozo et "aLI (|2007|) , we attempted to calibrate the 
various nuisance parameters in our model using mock 
catalogs created with the method of Wechsler et al. 2007, 
in preparation). We found that for some parameters, sys- 
tematic variation between simulations dominated over 
random errors. Specifically, while both the complete- 
ness and purity functions appeared to be stable and ro- 
bustly constrained, the parameters which characterizes 
the normalization of the probability matrix P(N {, s \Nt) 
and its scatter, B and B\ were seen to have large sys- 
tematic variations amongst the three mock catalogs con- 
sidered. These systematic variations, however, appeared 
to fall along a degeneracy band, as see n in the low e r pane l 
of Figure 5 of the companion paper IRozo et all lj2007h . 
and reproduced here in Appendix [XJ Consequently, we 
placed a weak prior on Bq and B\ along this degener- 
acy, and wide enough that it comfortably encompassed 
the 95% statistical region s of the parameter s in all three 
simulations considered in IRozo et al.l (|2007t ). Our prior 
is shown with solid lines in Figure [SJ The slope (3 of the 
mean relation (A {, s |Ay) appeared to be somewhat bet- 
ter constrained at roughly the ~ 5% — 10% level. In an 
effort to be conservative, and given the small number of 
simulations we had available, we have opted for placing a 
more generous 15% Gaussian prior on (3 centered on the 
simulation-calibrated value. Since both the completeness 
and purity functions appear to be robustly constrained 
in the simul ations, we a l so use the simulation-calibrated 
priors from IRozo et "all (|2007f ) for these quantities. Fi- 
nally, whereas we found the range of photometric rcdshift 

14 An example of such a distribution is Poisson statistics, in 
which case the proportionality constant is simply unity. 
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Fig. 2. — Comparison between the observed binned cluster counts 
in the maxBCG catalog (solid circles) and our maximum-likelihood 
model. Error bars represent Poisson uncertainties, and are shown 
for reference only as the counts in the various richness bins are 
correlated. The agreement between our maximum likelihood model 
and the observed number counts is excellent. 

parameters to also be systcmatics dominated, it was clear 
that the range of these uncertainties was small enough 
that photometric redshift uncertainties have a minimal 
impact in our results. Thus, we simply marginalize over 
the range of photometric rcdshift parameters observed in 
the simulations. 

Unfortunately, the above priors are simply not restric- 
tive enough for constraining cosmological parameters. In 
particular, the range of cosmological and HOD models is 
large enough that evaluation of the likelihood function 
over the entire degeneracy region becomes impossible. 
To overcome this difficulty, we therefore include three 
additional priors. The first is a CMB based prior on the 
matter density, £l m h 2 = 0.128 ± 0.01. It is important to 
note, however, that the CMB constraint on the matter 
density of the universe is independent of the amplitude 
of the power spectrum, as it depends only on the well 
known physics of the photon-baryon fluid of the early 
universe. Consequently, use of this prior in our cosmolog- 
ical analysis should not introduce a bias in our estimate 
for erg regardless of the dynamical nature of dark energy. 
In a similar spirit, we assume a supernova-based prior 
h = 0.73 ± 0.05, and finally, a generous theoretically- 
motivated prior on a, the slope of the HOD, which w e 
take to be a = 1.0 ±0.15 (see e.g. lKravtsov et al J 120041 ) . 
All of these priors are assumed to be Gaussian in log 
space (i.e. lognormal). 

To summarize, then, we have placed priors on the phys- 
ical parameters Q m h 2 , h, and a, the nuisance parameter 
/3, and the purity and completeness functions. The am- 
plitude and scatter of the A {, s — Nt relation are allowed 
to float essentially freely The remaining physical pa- 
rameters of interest are Mi, the mass scale of the HOD, 
and cTg , the amplitude of the power spectrum in cluster 
scales. 

4. RESULTS 

Constraints on our model parameters are obtained 
from the likelihood function described above via a Monte 
Carlo Markov Chain (MCMC) appro ach. The details o f 
our MCMC algorithm can be found in lRozo etldl (l2007h. 
Briefly, leaning heavily on the work by iDunklev et al.l 
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Fig. 3.— Confidence regions (68% and 95%) in the a 8 - Mi 
plane. The circle marks the median values of Mi and erg (8.4 ■ 
10 12 Mq, 0.92), while the diamond marks the corresponding aver- 
age values (9 ■ 10 12 Mq, 0.92). and Mi are the only two com- 
pletely free parameters in our analysis. 

(|2005f ). we construct MCMCs that arc optimal in step 
size, and we ensure we make enough evaluations to ro- 
bustly recover the 95% confidence likelihood contours in 
parameter space. The data itself is binned in nine log- 
arithmic bins in the richness range 100 > N f, s > 10, 
plus an additional high richness bin containing all clus- 
ters with 100 galaxies or more. The rcdshift range is 
[z m i n , z max ] = [0.1,0.3]. These bins are wide enough that 
even our least-populated bin contains 14 clusters, which 
is necessary given our likelihood function 15 . Finally, in 
order to make sure our results are robust, we ran two 
additional MCMCs, and checked that the recovered dis- 
tributions were consistent with each other. We found this 
to be the case. Consequently, we then proceeded to join 
all three chains into a single chain with 3T0 5 points. This 
number of evaluations is enough to recover the 95% confi- 
dence regions of the distribution with « 5% accuracy, or 
alternatively the 99% confidence regions with sa 10% ac- 
curacy. Finally, in order to test whether our results were 
robust to the number of clusters with most extreme rich- 
nesses, we also ran an MCMC where the richness range 
was limited to 100 > N a i, s > 11, which amounts to drop- 
ping the 14 richest clusters and the 2558 clusters with 
N i, s = 10. We found that this chain produced results 
consistent with those of our original analysis. 

A comparison between the observed number counts 
and the recovered maximum-likelihood model can be 
seen in Figure [3] We can see that our best fit model 
provides an excellent fit to the observed number counts. 
The corresponding 68% and 95% confidence regions in 
the a 8 — Mi plane are seen in Figure [3J The median 
and average values for each of these two parameters are 
(Mi,ct 8 ) = (8.4 • 1O 12 M , 0.92) and (Mi,a 8 ) = (9.0- 
10 12 M Q ,0.92) respectively. The marginalized distribu- 
tion for a& , shown in Figure[51 is roughly fit by a Gaussian 
distribution with eg = 0.92 ±0.10, and allows us to place 
an interesting lower limit on as- c 8 > 0.76 (95% CL) 
or a 8 > 0.68 (99% CL). The Mi distribution is roughly 

15 Specifically, the likelihood approximates the Poisson uncer- 
tainty in the abundance as Gaussian, so a large number of clusters 
per bin is necessary in our analysis 
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Fig. 4. — The posterior distribution in the a — f) plane. The filled 
contours represent the 68% and 95% confidence regions, which arc 
centered on the square; also shown for reference are the correspond- 
ing regions in our prior distribution. The circle marks the median 
values of each of the parameters (f3, a) = (1.15, 0.84), which are al- 
most identical to the average values. Note that the central value for 
our input priors (the square) lies outside the 95% confidence region 
of the a posteriori distribution, suggesting that either our calibra- 
tion is incorrect, or the slope of the HOD is substantially lower than 
unity. The solid line corresponds to the expected a/3 constant 
degeneracy between the two parameters. 

lognormal, though with some slight skewness. The cor- 
responding 1 — a parameters are ln(Mi/10 12 MQ) = 
2.1 ± 0.4, corresponding to Mi = 8.2^ ? x 1O 12 M . 

While our remaining model parameters were all con- 
strained with priors, it is nevertheless interesting to look 
at their a posteriori distributions. In particular, we find 
that there appears to be some tension between our prior 
on the slope (3 of the N b s — Nt relation, and our prior 
on a, the slope of the HOD. This is illustrated in Fig- 
ure 31 where we show the posterior distribution in the 
a — [3 plane. The circle marks the median values of the 
marginalized distributions (/3, a) = (1.15,0.84) whereas 
the square marks the central values of our priors, a = 1.0 
and (3 = 1.18. The average values of the parameters are 
very close to the median values. The degeneracy direc- 
tion is roughly a[3 = constant, as expected based on the 
fact that N obs ~ N% ~ m afi . While the width of our 
priors is large enough that we do not feel this discrep- 
ancy biases our results, it is clear from Figure 2] that 
either the slope of the N b s — Nt relation is shallower 
than seen in the simulations, or the slope of the HOD is 
markedly different from unity. Information from a new 
suite of simulations for currently ongoing work indicates 
that it is probably the former: our prior for (3 seems to 
be on the high edge of what is possible. 

While it is impossible for us to fully determine whether 
the HOD prior or our simulation-calibrated prior is more 
incorrect without additional information, we can try to 
better understand how each of these two possibilities 
would affect our results. To do so, we have run two 
additional MCMCs, one with a tight, 5% prior on a and 
no prior on /?, and one with the converse priors. Due to 
the strong a(3 degeneracy, in either case we found that 
the maximum likelihood model resulted in an excellent 
fit to the data. 

Our results are summarized in Figure [5j Briefly, we 
find that the tight (3 prior favors a low a [a = 0.76 ± 
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Fig. 5. — The solid histogram shows the marginalized erg dis- 
tribution recovered from our MCMC. A Gaussian fit to the distri- 
bution results in erg = 0.92 ± 0.10. From the above distribution 
we can also directly obtain a lower limit for erg, erg > 0.76 at 95% 
confidence, or erg > 0.68 at 99% confidence. Also shown are the 
marginalized distributions recovered when one places a tight, 5% 
prior on the slope of the HOD (dotted line) and no prior on the 
selection function, and for the converse case (dashed line). 

0.05) and a very high a$ (as = 1.05lg Jj), whereas the 
tight HOD prior a = 1.0 ± 0.05 favors a lower as value. 
a% = 0.92 ± 0.11. The implausibly large value for as 
and the further lowering of a for the (3 prior suggests 
that our central value for (3 may be too high, which is 
consistent with more recent preliminary determinations 
of the maxBCG selection function in simulations of our 
ongoing work. 

5. DISCUSSION 

5.1. Comparison With Other Work 

We have performed a detailed statistical analysis of 
the ob served cluster counts i n the maxBCG cluster cat- 
alog of iKoester etHI (|2007bf l. We find cr 8 = 0.92 ± 0.10 
(1-er), and place an interesting lower limit rjg > 0.75 
(95% CL) or as > 0.66 (99% CL). All of our results are 
marginalized over all major systematic effects, though 
they are also subject to additional cosmological priors 
n m h 2 = 0.128 ± 0.01 and h = 0.73 ± 0.05. Finally, our 
analysis was restricted to flat cosmologies with massless 
neutrinos and scale-invariant primordial power spectra. 

Our value for as is somewhat larger — but con- 
sistent with — the typical value obtained from clus- 
ter abundance studies using X-rays (as y 0.7, 
see e.g. iBorgani et all 120011: iPierpaoli etaTl 1 2001 



Ikebe et al.ll2002tlAllen et al.ll2003tlPierDaoli et afl l2003 



Schuecker et al. | 120031 : iHenrvl 120041 . and references 
therein), as determinations from optical cluster samples 
are more rare, but here again we find our results to be in 
good agreeme nt with previous work. For instance , both 
iBerlind et ail (|2006l ) and | van den Bosch et alj (12006) 
suggest as ~ 0.8, though no errors on t his es timate are 
made. On the other hand, iRines et al.l (|2006f ) analyze a 
local cluster sample selected from the SP SS and report 



a 8 = 0.92 



+0.24 
-0.19' 



similar to the result from iBahcall et alJ 
(l2003h and that of lBahcall fc Hodel (|200l . 

The work that is per haps most direct l y com parable in 
spirit to ours is that of Idadders et all (|2006l ). who find 



lyze cluster abundances as observed in the Red-Sequence 
Survey (RCS). As in our work, their value for as is 
marginalized over uncertainties in both the mass-richness 
relation and the hubblc parameter h. This is particu- 
larly important as it is known that these are the param- 
eters that are m ost highly degenerate with as (see e.g. 
iRozo et all l2004h . yet they are often held fixed in clus- 
ter abundance studies. Despite the fa ct that the central 
value for a s in Idadders et al.l (|2006r i is below our 99% 
confidence lower limit as > 0.68, given the large error 
bars in both of our determinations it is clear that our re- 
sults are, in fact, quite consist ent with each other. M ore- 
over, the likelihood function in lGladders et alJ (|2006f ) has 
a large tail that extends out to erg ps 1.2, implying the 
overlap between the two analysis is even larger than what 
the 68% confidence error bars might suggest. 

What is clear at this time is that a precise determina- 
tion of as still eludes us. This is unfortunate, as our large 
errors for as imply that almost any reasonable value that 
wc could find would automatically be c onsistent with 
the K CDM interpolated value found by ISpergel et al.l 



(Ell), 



tT8 



0.74 



+0.05 
-0.06- 



It is evident that an improve- 



rs = 0.67 



+0.18 
-0.13 



using a self-calibration technique to ana- 



mcnt of at least a factor of two in the error bars is needed 
before as constraints from cluster abundances become a 
potential probe of the evolution of dark energy Fortu- 
nately, such an improvement may be possible in the near 
future. For instance, we are now in the process of measur- 
ing weak lensing masses for cl usters of fixed richness N b s 
using the method presented in lJohnston et all (|2004f ) . In 
principle, such a measurement represents a direct mea- 
surement of the richness-mass relation, and hence can 
place direct constraints on the mass scale Mi and the 
product a/3, the effective slope of the mass-richness rela- 
tion. 

An additional interesting result at this stage is the low 
value for a, the slope of the HOD, that we recovered 
with our analysis. Whereas we find a = 0.83 ± 0.06, 
there is a significant body of evidence that points to- 
wards larger values of a. From the observational side, 
iKochanek et alJ (|2003l ) obtain a slope a ~ 1 or a lit- 
tle steeper on the basis of a 2MASS cluster sample 
obtained from a matched filter algorithm, though we 
note that a recent new analysis based on stacked X- 
ray emission fro m these cluster s suggests a lower value 
a = 87 ±0.05 (jDai et al.ll2006l) . Likewise. [Tinker et all 
(|2006f ) find a slope of unity based on the galaxy angular 
correlation fun ction and the distrib ution of voids in the 
SDSS. Finally, IZehavi et all (|2005h find that the slope 
a is close to unity, but steadily increases with increas- 
ing luminosity. They also note, however, that this re- 
sult is dependent on the detailed form of the HOD pa- 
rameterizati on. In particular, an alternate parameter- 
ization from iKravtsov et alJ (|2004[ ) with slope of unity 
on the high mass end is seen to be con sistent with the 
data a s well. From the theoretical side, IKravtsov et alJ 
(|2004| ) showed that the slope a characterizing the HOD 
of subhalos is very close to unity. Moreover, numeri- 
cal simulations of galaxy formation suggest there is an 
excellent corr espondence betw e en dark matter subhalos 
and galaxies dGao et all 120041: iNagai fc Kravtsovl [20051 : 
IWeinberg et all l2006f ). This is further evidenced by 
the fact that simple models used to assign luminosities 
to galaxies provide excellent matches to the luminosity 
dependent galaxy-galaxy and galaxy-mass correlations 
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functions (iTasitsiomi et al.l l2004t Mandcl baum et all 
[20MlConrov et all l200ftl Vale & Ostrikerll200ffl . as well 
as the galaxy three-point function (Marin et al. 2007, in 
preparation) . 

There are, however, some lines of evidence for a lower 
value for a. In particular, some numerical simula- 
tions suggest a valu e a tts 8 on the high mass end 
(|Berlind et al.l l2003f) . Also, IScoccimarro et al.l (|200lD 
cited a low value of the slope a = 0.8 on the basis of the 
galaxy three-point function, though we note they did not 
differentiate between central and satellite galaxies, which 
would tend to l o wer t he recovered slope for the HOD. 
lAbazaiian et al.l (|2005f ) also found that the best fit value 
for a was a = 0.83jlo' 2 3 on the basis of a joint analysis 
of the galaxy angular correlation function and the CMB, 
though note the rather large error bars. Finally , analy - 
sis of the 2dFGRS cluster catalog of lYang et all (j2005bl ) 
also suggests that if a$ w 0.9 then the number of galax- 
ies in massive halos is low relative t o standard models , 
consistent with a lower value for a (|Yang et al.ll2005ah . 
Nevertheless, they emphasize their data is equally well 
fit by a low as model with a standard galaxy population, 
and in their most recent analysis they argue that a mod- 
erately low da (roughly 0.8 < as 5, 0.9) is likely to be 
correct (jvan den Bosch et al.ll2006( ). 

Overall, the majority of the evidence available suggests 
that the slope of the HOD is typically closer to unity 
than what we have found. Of course, in general one ex- 
pects that in detail, the slope of the HOD will depend 
on the particular criteria for galaxy selection in comput- 
ing the HOD. At this point, about all we can say is that 
the maxBCG cluster sample provides some marginal ev- 
idence for a < 1, though this is only about a 2tr effect, 
and, moreover, the deviation from unity might be slightly 
overestimated if indeed the selection function in the sim- 
ulations and in the real data are different. 

5.2. Additional Sources of Systematics 

An absolutely fundamental assumption about our sta- 
tistical method is that the cluster selection function is an 
inherent property of the cluster-finding algorithm used. 
More specifically, we have assumed that given a halo of 
richness Nt, the cluster-finding algorithm has a proba- 
bility P(N i,s\Nt) of detecting such a halo as a cluster 
of richness N a b s ■ If this probability is not a property of 
the cluster-finding algorithm (ie, if it is a strong function 
of cosmology), then our analysis needs to be generalized, 
and the cosmological dependence of the probability ma- 
trix needs to be calibrated. Consequently, the extent 
to which the above probability is robust to moderate 
changes in cosmology is an inherent systematic of our 
method, and clearly warrants further investigation. We 
emphasize, however, that this is true of any characteriza- 
tion of a cluster selection function obtained through the 
use of simulations. That is to say, it is not guaranteed 
a priori that two different yet realistic simulations will 
result in the same cluster selection function. 

In addition to the above systematic, by far our most 
important uncertainties are due to possible selection ef- 
fects not included within our model. For instance, we 
expect each of the components of our model — the com- 
pleteness and purity functions and the signal matrix — 
to have some redshift and richness dependence. In this 
work, we have simply ignored this possibility, though we 



note that due to the small range of redshifts probed, we 
do not expect this systematic to be particularly signif- 
icant. Moreover, we have explicitly demonstrated that 
ignoring such evolution in the simulation s does not bias 
our results in any way (R ozo "etld1l200l . 

In addition to these selection function systematics, 
there are additional sources of error which we have not 
included. For instance, in this work we have ignored pos- 
sible evolution in the HOD. Given the relatively narrow 
redshift range considered, and the fact that there appears 
to be little evolution in the way galaxi es populate halo s 
between redshifts z = and z = 0.8 (|Yan et al.ll2003h . 
we do not expect the no evolution assumption to be a 
limiting factor in our analysis. 

A more theoretical systematic which we have not con- 
sidered has to do with the current uncertainty in the 
predicted halo mass function. In particular, while the 
halo mass function appe ars to be universal w ith about a 
« 20% margin of error ([Jenkins et al.ll200ll ). additional 
work is required to test whether the halo mass function 
is indeed universal to higher accuracy — and if not, to 
characterize any intrinsic cosmological dependences. In 
this work, we have ignored these complications and have 
made no attempt to marginalize over the correspond- 
ing uncertainties (as is customary in the literature). We 
have, however, explicitly checked that changing the pa- 
rameterization of th e halo ma ss function from t hat o f 
I Jenkins et all (I2001D to that of IShet h fc Tormenl (|2002f ) 
or that of I Warren et al.l (|2005l ) changes the expectation 
values of our parameters by much less than our quoted 
la uncertainty. 

Yet another possible systematic of theoretical origin 
has to do with our assumptions about the HOD. In par- 
ticular, if the HOD has any curvature over the mass 
range probed, then our model will necessarily result in 
biased parameter estimation. Since at this time we are 
only probing roughly one and a half decades in mass, 
we do not expect this to be a significant problem. Nev- 
ertheless, once the selection function for the maxBCG 
cluster-finding algorithm is better understood, it would 
be interesting to investigate to what extent the data can 
constrain deviations from linearity in the mean relation 
between halo mass and Nt- 

5.3. Future Work and Improvements 

Clearly, one of the most important problems to work 
on at this time is improving our understanding of the 
maxBCG selection function. Note that this work must 
involve investigating a large range of cosmological pa- 
rameters so that any cosmology dependences inherent 
to the cluster-finding algorithm can be adequately cali- 
brated. In addition, it is possible the simulations them- 
selves need to be refined to produce more realistic skies so 
that differences in the cluster selection function between 
the simulations and the real sky are minimized. 

Along this same line of reasoning, an important possi- 
bility that must be considered is to ask what the most 
useful richness definitions are, and its related question, 
what the best way of matching halos to clusters is. While 
we have do ne s ome p reliminary work in this direction 
iRozo et "all (see I2007T ). we note that membership-based 
matching algorithms obviously depend on what we mean 
by a halo member. Thus, if the definition for Nt changes, 
the matching of halos to clusters, i.e. the selection func- 
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tion, changes. Indeed, preliminary work suggests that 
uncertainties in the selection function of the cluster- 
finding algorithm are significantly reduced if only ridge- 
line galaxies are considered as halo members. Note that 
from a theoretical perspective, such a definition for Nt 
would be perfectly fine, since galaxy formation models 
suggest that early ty pe galaxies in clus ters also obey a 
simple Poisson HOD (jZheng et alj |2005). Of course, the 
important thing is not the fact the HOD is Poisson, but 
that we have a theoretically well-motivated reason for 
choosing a particular form for the scatter between halo 
mass and richness. 

Finally, it is evident that perhaps the single most im- 
portant step at this time will be to include additional 
data that will soon be available for the maxBCG sam- 
ple. Specifically, we are now in the process of computing 
ensemble-averaged X-ray, dynamical, and weak lensing 
masses for the SDSS maxBCG cluster sample. These 
new data sets will directly constrain the richness-mass 
relation of the clusters, and their inclusion will provide 
an entirely new tool that will further improve the cosmo- 
logical constraints derived here. 

In light of the above discussion, it is clear that there 
remains an enormous amount of work to be done to fully 
realize the potential of optical cluster catalogs. Never- 
theless, we emphasize that, even in this first, roughest 
attempt to recover science from the SDSS maxBCG clus- 
ter sample, we have been able to derive robust cosmo- 
logical constraints that are competitive with other ap- 
proaches, squarely placing optical cluster science in the 
general toolkit of the precision cosmology effort. 
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APPENDIX 

DEFINITION OF THE NUISANCE PARAMETERS 

Here we summarize the expressions that define th e nuisance param eters in our model. For a discussion of where 
these expressions come from, we refer the reader to iRozo et al.l (|2007j ). The probability distribution P(N obs \NT) is 
characterized as follows: let (j,ml(Nt) be the most likely value for N obs given Nt- (J-ml is parameterized as a power 
law 

Vml(Nt) = 20exp(B )(AW20) /3 , (Al) 
which defines Bq and (3. The mean value of N obs at fixed Nt is found to be given by 

{N obs \N T ) = 20cxp(B + 0.14)(A^ T /20)(' 3 - - 12 ) (A2) 

and the variance is given by 

Vzr(N obs \N T ) = cxp(-3B + Bi)/x(iV T ) (A3) 

which defines B\. The confidence regions in the Bq — B\ plane for each of the three simulations considered in lRozo et al.l 
(2007) are seen in figure [6] The solid lines define a band over which the parameters Bq and B\ are allowed to vary. 

The completeness function is observed to be richness independent, and is thus parameterized in terms of a single 
parameter c such that c(Nt) = c. The purity function is characterized with two parameters po and p± such that 

p(N obs ) = cxp(~x(N obs ) 2 ) (A4) 

where 

'^^^(^T)- 1 )- (A5) 

Finally, the photometric rcdshift distribution p(z c \zh) where z c is the observed photometric rcdshift of a cluster and 
Zh is the true spectroscopic redshift of the parent halo is taken to be of the form 

p{z c \zh) = —Pb{b\z h ) (A6) 

Zh 

where b = z c /z b is the photometric rcdshift bias. p b is found to be Gaussian, and the nuisance parameters (b z ) and 
a z correspond to the mean and standard deviations of said Gaussian. 
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FlG. 6. — 95% confidence regions for the Bo and B± in each of the three mock catalogs analyzed in Rozo et al.l (|2007h . The dashed 
and solid contours are for Mocks A and B respectively. The shaded contours are 68% and 95% confidence regions in Mock C. The small 
filled circles mark the best-fit parameters from the mock catalogs, and were used to generate the Monte-Carlo realizations from which the 
confidence regions are derived. The solid lines mark the band over which the parameters Bq and B\ arc allowed to vary. 



